      Program Stat
      !*****************************************************************
      include 'accur.for'
      include 'statd.for'
      Dimension r(0:IMax),alpha(0:IMax),theta(0:JMax)
      Common /phi/phi(0:IMax,0:JMax)
     +       /delta/delta
     +       /pi/pi
     +       /IB/IB/JB/JB/I0/I0
     +       /rho0/rho0/phiu/phiu
     +       /gamma_grav/gamma_grav
      Logical LogIter,LogInit,LogRapid

      pi=acos(-one)
      ! === constants ==>
      LogRapid=.true.
      JB=0 ! spheroidal configurations
      delta=1d0 ! grid parameter
      rho0=1d0
      gamma_grav=one
      ! <== constants ===
      open(1,file='output',access='append')
      Call timer(ITime0) ! PC
          write(1,*) '    rho0   ',
     +               'r(IB)**delta',
     +               '     CM    ',
     +               '     CJ    ',
     +               '  T/abs(W) ',
     +               '   gamma_  ',
     +               '     -W    ',
     +               '     VT     '
      close(1)
      dr=one/IA
      do i=0,IA
          r(i)=dr*i
      enddo
      dalpha=one/r(IA-1)-one/r(IA)
      do l=0,IMax-IA+1
          alpha(l)=dalpha*l
      enddo
      do i=IA+1,IMax-1
          l=IMax-i
          r(i)=one/alpha(l)
      enddo
      dtheta=half*pi/JMax
      do j=0,JMax
	  theta(j)=dtheta*j
      enddo
c     === initial phi ==>
          do j=0,JMax
              do i=0,IMax
                  if(i.le.IA) then
                      phi(i,j)=float(i-IA-1)/float(IA)
                  else
                      l=IMax-i
                      phi(i,j)=phi(IA,j)*degree(alpha(l),delta)
                  endif
              enddo
          enddo
c     <== initial phi ===
      if(JB.eq.0) then ! spheroidal configurations
          I0=0
      else ! toroidal configurations
          I0=half*(IB+IA)
      endif

      IB=IA
      do while(IB.ne.0) ! ===== IB cycle ====>

          LogInit=.true.
*         write(*,*) 'IB',IB
          C=( phi(IA,JMax)+Fpsi(one)
     +       -phi(I0,JMax)*EOS(zero,3)/EOS(rho0,3)
     +      )
     +      /(one-EOS(zero,3)/EOS(rho0,3))
          phiu=EOS(rho0,3)/(C-phi(I0,JMax))
	  n_max=20
	  LogIter=.true.
	  n=1
          do while(LogIter.and.(n.le.n_max))
	      C0=C
              Call new_phi(r,alpha,theta,dtheta,
     +                     C,phiu,rho0,LogInit,LogRapid)
	      if(dabs((C0-C)/C).lt.1d-7) LogIter=.false.
*              write(*,*) n,C,phi(I0,JMax)
	      n=n+1
          enddo
1         Format(i5,4(1pe10.2),1pe28.20)
          Call Integrals(CM,V,CJ,T,W,Eth,Pip,phi,r,theta,C,phiu,
     +                   rho0,omegaA)
	  VT=abs(2d0*T/W+one+3d0*Pip/W)
          rhomax=zero
          do i=0,IA
              do j=0,JMax
		  rhomax=max(rhomax,
     +                     EOS(phiu
     +                         *( C
     +                           -phi(i,j)
     +                           -Fpsi(degree(r(i),
     +                                        delta)*sin(theta(j)))),
     +                          1
     +                     )/rho0)
              enddo
          enddo
          open(1,file='output',access='append')
          write(1,2) rho0,r(IB)**delta,CM,
     +               CJ,
     +               T/abs(W),Pip/Eth+one,
     +               -W,
     +               VT
2         Format(10(1pe11.3))
          close(1)
          write(*,*) 'IB, VT', IB, VT
          Call Sigma(C)
          if(IB.ge.0.1d0*IA) then
              IB=IB-1
          else
              IB=0
          endif
      enddo
      Call Timer(ITime)
      open(1,file='output',access='append')
      write(1,*) '* time *',(ITime-ITime0)/1d2,'sec'
      close(1)
      stop
      end

      double precision Function EOS(x,n)
      !*****************************************************************
      include 'accur.for'
      Parameter(CN=1.5d0)
      Parameter(CK=1d0)
      Parameter(gamma=one+one/CN)
      if(n.eq.1) then ! rho(H)
          if(x.gt.zero) then
              EOS=( x*(gamma-one)/(gamma*CK)
     +            )**(one/(gamma-one))
          else
              EOS=zero
          endif
      endif
      if(n.eq.2) then ! drho_dH(H)
          if(x.gt.zero) then
                  EOS=(x*(gamma-one)/(gamma*CK)
     +                )**(one/(gamma-one)-one)
     +                /(CK*gamma)
          else
              EOS=zero
          endif
      endif
      if(n.eq.3) then ! H(rho)
          if(x.gt.zero) then
                  EOS= CK*gamma
     +                 *( x**(gamma-one)
     +                  )
     +                 /(gamma-one)
          else
              EOS=zero
          endif
      endif
      if(n.eq.4) then ! p(H)
          if(x.gt.zero) then
                  rho_=( x*(gamma-one)/(gamma*CK)
     +                 )**(one/(gamma-one))
              EOS=CK*rho_**gamma
          else
              EOS=zero
          endif
      endif
      if(n.eq.5) then ! E(rho)
          if(x.gt.zero) then
                  EOS= CK
     +                 *( x**(gamma-one)
     +                  )
     +                 /(gamma-one)
              EOS=EOS*x
          else
              EOS=zero
          endif
      endif
      return
      end
